function osl_test_script_sensorspace(testdatadir,testoutputdir, S)

testoutputdir=[testoutputdir '_sensorspace'];
runcmd(['rm -rf ' testoutputdir]);
mkdir(testoutputdir);

testplotsdir=[testoutputdir '/plots'];
mkdir(testplotsdir);

if(1),
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% osl_example_sensorspace_oat

printprefix='sensorspace_oat';
printindex=1;

datadir=[testdatadir '/faces_subject1_data']; % directory where the data is

clear spm_files_continuous spm_files_epoched;
% set up a list of SPM MEEG object file names (we only have one here)
spm_files_continuous{1}=[datadir '/spm8_meg1.mat'];
spm_files_epoched{1}=[datadir '/espm8_meg1.mat'];

oat=[];
oat.source_recon.D_continuous=spm_files_continuous;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
oat.source_recon.D_epoched=spm_files_epoched; % this is passed in so that the bad trials and bad channels can be read out
oat.source_recon.freq_range=[]; % frequency range in Hz
oat.source_recon.time_range=[-0.2 0.4];
oat.source_recon.method='none';
oat.source_recon.dirname=[testoutputdir '/' printprefix];

Xsummary={};Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=Xsummary;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
oat.first_level.contrast{4}=[0 0 -1 1]'; 
oat.first_level.contrast{5}=[0 -1 0 1]'; 
oat.first_level.contrast_name{1}='motorbikes';
oat.first_level.contrast_name{2}='faces';
oat.first_level.contrast_name{3}='faces-motorbikes';
oat.first_level.contrast_name{4}='fear-happy';
oat.first_level.contrast_name{5}='fear-neutral';

oat.first_level.cope_type='acope';

oat = osl_check_oat(oat);

oat.to_do=[1 1 0 0];
oat = osl_run_oat(oat);

% visualise using Fieldtrip
S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.modality='MEGPLANAR';
S2.first_level_contrast=[1 2 3];
S2.cfg.interactive = 'no';

% calculate t-stat using contrast of absolute value of parameter estimates
osl_stats_multiplotER(S2);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% osl_example_sensorspace_oat_tf

printprefix='sensorspace_oat_tf';
printindex=1;

%%%
% DO SENSOR SPACE MULTI-BAND TIME-FREQ ANALYSIS USING OAT
 
datadir=[testdatadir '/faces_subject1_data']; % directory where the data is

clear spm_files spm_files_epoched;
% set up a list of SPM MEEG object file names (we only have one here)
spm_files{1}=[datadir '/spm8_meg1.mat'];
spm_files_epoched{1}=[datadir '/espm8_meg1.mat'];

oat=[];
oat.source_recon.D_continuous=spm_files;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
oat.source_recon.D_epoched=spm_files_epoched; % this is passed in so that the bad trials and bad channels can be read out
oat.source_recon.freq_range=[2 30]; % frequency range in Hz
oat.source_recon.time_range=[-0.2 0.4];
oat.source_recon.method='none';
oat.source_recon.dirname=[testoutputdir '/' printprefix];

oat.first_level.tf_method='hilbert'; % can be morlet or hilbert
oat.first_level.tf_num_freqs=8; % we are keeping this unusally low in the practical for the sake of speed
oat.first_level.tf_hilbert_freq_res=4;
oat.first_level.bc=[1 1 0];

Xsummary={};Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=Xsummary;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes

oat = osl_check_oat(oat);

oat.to_do=[1 1 0 0];
oat = osl_run_oat(oat);

% load GLM result
stats_hilb=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% visualise using Fieldtrip
S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.modality='MEGPLANAR';
S2.first_level_contrast=[3];
S2.cfg.colorbar='yes';
S2.cfg.interactive = 'no';

% calculate t-stat using contrast of absolute value of parameter estimates
[cfg, data]=osl_stats_multiplotTFR(S2);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%
% DO SENSOR SPACE SINGLE-BAND TIME-FREQ ANALYSIS USING OAT

printprefix='sensorspace_oat_singtf';
printindex=1;

datadir=[testdatadir '/faces_subject1_data']; % directory where the data is

clear spm_files spm_files_epoched;
% set up a list of SPM MEEG object file names (we only have one here)
spm_files{1}=[datadir '/spm8_meg1.mat'];
spm_files_epoched{1}=[datadir '/espm8_meg1.mat'];

oat=[];
oat.source_recon.D_continuous=spm_files;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
oat.source_recon.D_epoched=spm_files_epoched; % this is passed in so that the bad trials and bad channels can be read out
oat.source_recon.freq_range=[5 20]; % frequency range in Hz
oat.source_recon.time_range=[-0.1 0.3];
oat.source_recon.method='none';
oat.source_recon.dirname=[testoutputdir '/' printprefix];

oat.first_level.tf_method='hilbert'; % can be morlet or hilbert
oat.first_level.tf_num_freqs=1; % we are keeping this unusally low in the practical for the sake of speed
oat.first_level.tf_hilbert_freq_res=diff(oat.source_recon.freq_range);

Xsummary={};Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=Xsummary;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes

oat = osl_check_oat(oat);

oat.to_do=[1 1 0 0];
oat = osl_run_oat(oat);

% load GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.modality='MEGPLANAR'; 
S2.first_level_contrast=3;
S2.cfg.interactive = 'no';

% calculate t-stat using contrast of absolute value of parameter estimates
[cfg, data]=osl_stats_multiplotER(S2);

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% osl_example_sensorspace_continuous_oat

printprefix='sensorspace_continuous_oat';
printindex=1;

datadir=[testdatadir '/ctf_fingertap_subject1_data']; % directory where the data is

% Set up the list of subjects and their structural scans for the analysi 
% Currently there is only 1 subject.
clear spm_files;

% set up a list of SPM MEEG object file names (we only have one here)
spm_files={[datadir '/dsubject1.mat']};
 
oat=[];
oat.source_recon.D_continuous=spm_files;
oat.source_recon.conditions={'Undefined'};
oat.source_recon.freq_range=[13 30]; % frequency range in Hz
oat.source_recon.time_range=[300,32*30];
%oat.source_recon.time_range=[300,14*30];
oat.source_recon.method='none';
oat.source_recon.modalities={'MEGGRAD'};
oat.source_recon.work_in_pca_subspace=1;
oat.source_recon.pca_dim=270;

oat.source_recon.dirname=[testoutputdir '/' printprefix];

oat = osl_check_oat(oat);

oat.to_do=[1 0 0 0];

oat = osl_run_oat(oat);

%%%
% Establish regressor for continuous time GLM. 

D=spm_eeg_load(spm_files{1});

% This should be setup to correspond to the same time window
% as the full time window for D

x=zeros(length(D.time),5);

block_length=30; %s
block_order=[5 5 5 5 5 5 5 5 5 5 4 3 2 1 2 3 1 4 3 4 1 3 2 1 4 4 2 1 3 3 4 1 4 3 1 2 1 2 3 4 3 4 1 2 3 4 1 2];

% [Left, Right, Rest, Both, Rest_at_start]
% [  1     2      4     8     16 ]
% figure;plot(D.time,squeeze(D(1,:,:)))
% emacs ~/vols_data/From_Nottingham_with_Love/JRH_MotorCon_20100429_01_FORMARK.ds/MarkerFile.mrk 

tres=1/(D.fsample);
tim=1;
for tt=1:length(block_order),    
    x(tim:tim+block_length/tres-1,block_order(tt))=1;
    tim=tim+block_length/tres;
end;

figure;plot(D.time,x);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%%
% Run time-wise glm to do regression against known finger tapping regressors

oat=osl_load_oat(oat.source_recon.dirname);

oat.first_level.tf_method='hilbert';
oat.first_level.tf_downsample_factor=10; 
oat.first_level.time_moving_av_win_size=1;
oat.first_level.design_matrix=x';
oat.first_level.contrast{1}=[-1 0 1 0 0]'; % rest-left
oat.first_level.contrast{2}=[0 -1 1 0 0]'; % rest-right
oat.first_level.contrast{3}=[0  0 1 -1 0]'; % rest-both
oat.first_level.freq_range=[13 30];
oat.first_level.tf_hilbert_freq_res=diff(oat.first_level.freq_range);
oat.first_level.doGLM=1;

oat.first_level.name=['subj1_first_level_ft'];

oat.to_do=[0 1 0 0];

oat = osl_run_oat(oat);

% visualise using Fieldtrip
S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.modality='MEG';
S2.first_level_contrast=1;
S2.cfg.colorbar='yes';
S2.cfg.interactive='yes';
S2.cfg.zlim=[0 40];
%S2.cfg.interactive = 'no';

% calculate t-stat using contrast of absolute value of parameter estimates
osl_stats_multiplotER(S2);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

% try taking a look at the other contrasts and look at the lateralisation

%%%
% Look at beta power time course at a sensor over motor cortex

stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
% chan_label='MRF54';chan_ind=find(strcmp(stats.D_sensor_data.chanlabels,chan_label));
[a chan_ind]=max(squeeze(stats.cope(:,1,S2.first_level_contrast)))

% re-run oat without doing GLM - this will output the time series used to
% fit the GLM to, i.e. the beta power (Hilbert envelope) time courses
oat.to_do=[0 1 0 0];
oat.first_level.doGLM=0;
oat = osl_run_oat(oat);

% do plot
stats2=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
figure;plot(stats2.glm_input_times,normalise(squeeze(stats2.glm_input_data(chan_ind,:))));
% compare to design matrix:
time_ind=intersect(find(D.time>=oat.source_recon.time_range(1)),find(D.time<=oat.source_recon.time_range(2)));
ho;plot(D.time(time_ind),x(time_ind,1),'r','LineWidth',2);
plot(D.time(time_ind),x(time_ind,2),'g','LineWidth',2);
plot(D.time(time_ind),x(time_ind,4),'k','LineWidth',2);
legend('data','left','right','both');
plot4paper('time(secs)','beta power');

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

disp('*************************************************');
disp('Finished osl_test_scrip_sensorspace test');
disp('*************************************************');